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Progress  Report 


1.1  Introduction 


X-ray  diffraction  contrast  tomography  (DCT)  is  a  nondestructive  characterization  method  providing 
access  to  the  three-dimensional  grain  microstructures  in  a  wide  range  of  polycrystalline  materials.  DCT 
is  a  truly  three-dimensional  tomographic  imaging  approach,  sharing  a  common  experimental  setup  with 
conventional  X-ray  microtomography.  After  interaction  with  the  material,  both  the  transmitted  and 
diffracted  beams  are  recorded  on  a  high-resolution  X-ray  imaging  detector  positioned  close  to  the  sample. 
The  three-dimensional  distribution  of  the  X-ray  attenuation  coefficient,  and  the  three-dimensional  shape, 
grain  average  orientation  and  elastic  strain  tensor  of  all  grains  in  the  illuminated  sample  volume,  are 
determined  from  analysis  of  the  transmitted  and  diffracted  intensities,  respectively 

DCT  analysis  of  polycrystalline  materials  is  currently  a  hot  topic  at  synchrotron  radiation  facilities 
around  the  globe.  The  monochromatic  nature  of  a  synchrotron  beam,  its  parallel  beam  geometry, 
and  its  high  spatial  resolution  make  it  an  ideal  imaging  modality  to  measure  three-dimensional  grain 
microstructures  in  a  wide  range  of  polycrystalline  materials.  The  combination  of  three-dimensional  X-ray 
imaging  and  diffraction  techniques  on  the  same  instrument  are  particularly  interesting  for  the  study  of 
deformation  and  damage  mechanisms  in  structural  materials  [1,  2,  3],  but  current  applications  cover  a 
much  wider  spectrum,  including  phase  transformations  and  grain  coarsening  processes. 

Fig.  1.1  showcases  the  combined  use  of  (synchrotron)  X-ray  phase  contrast  and  X-ray  diffraction  con¬ 
trast  tomography  for  characterization  of  fatigue  crack  propagation  in  metal  alloys.  The  3D  dimensional 
grain  microstructure  of  a  bcc  Ti  alloy  sample  has  been  mapped  out  with  DCT  prior  to  an  in-situ  fatigue 
experiment  carried  out  on  the  same  specimen.  The  combined  characterization  enables  characterization 
and  interpretation  of  crack  orientation  and  crack  growth  rate  as  a  function  of  local  crystallographic 
orientation  [2]. 

However,  given  the  high  operating  cost  and  their  limited  availability,  synchrotron  facilities  often 
do  not  provide  a  feasible  solution.  It  would,  therefore,  be  very  beneficial  if  DCT  imaging  were  to  be 
possible  on  conventional  table-top  /iCT  systems  that  are  relatively  cheap  and  more  widely  available.  The 
downside  of  such  a  lab  source  CT  system  is  its  cone  beam  geometry  and  the  polychromatic  nature  of  its 
X-ray  source.  This  results  in  complicated  software  simulations  of  diffraction  projections  (and  thus  also 
their  reconstruction),  which  are  difficult  to  implement  efficiently  on  modern  computational  infrastructure 
such  as  GPU’s. 

In  this  work,  an  accurate  and  efficient  projection  model  is  introduced  for  the  case  of  polychromatic 
cone  beam  DCT.  In  Section  1.2,  the  forward  and  backprojection  model  of  this  projector  is  mathematically 
introduced  and  its  efficient  implementation  is  discussed.  Subsequently,  in  Section  1.3,  these  projectors 
are  used  to  reconstruct  a  grain  of  an  Aluminium  Copper  microstructure,  scanned  in  a  /r CT  system. 
Section  1.4,  finally,  concludes  this  work  and  provides  an  outline  into  future  research  directions  that  could 
greatly  benefit  the  application  of  DCT  in  lab  source  systems. 
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1.2.  CONE  BEAM  DCT  PROJECTOR 


Figure  1.1:  (a-c)  3D  rendition  of  a  fatigue  crack  in  a  Ti  alloy  reconstructed  from  phase  contrast  tomography  (PCT) 
after  46,  61  and  75  *103  loading  cycles.  The  color  corresponds  to  the  vertical  crack  position,  (d)  3D  rendition  of 
the  combined  data  set  showing  the  crack  and  the  surrounding  3D  grain  microstrucutre  as  determined  from  DCT.  (e) 
Surface  mesh  representing  the  fracture  surface,  colour  coded  with  respect  to  its  crystallographic  orientation.  Grain 
boundaries  are  labelled  in  white.  (f,g)  Crack  propagation  within  the  grain  labelled  in  (e).  The  activated  (110)  slip 
plane  and  the  two  (111)  slip  directions  are  highlighted  in  gray  and  black,  respectively. 


1.2  Cone  beam  dct  projector 

In  this  section,  a  projector  is  introduced  that  simulates  the  diffraction  images  of  grains  in  a  lab  source 
cone  beam  /rCT-system.  We  assume  that  the  diffraction  spots  have  already  been  indexed  and  that  for 
each  of  the  grains  the  center  of  mass  and  crystallographic  orientation  has  been  determined.  For  details 
on  the  indexing  procedure  for  polychromatic  DCT  data  we  refer  to  the  article  by  [4].  Note  that  the 
number  of  projection  images  per  grain  depends  on  material  and  experiment  parameters  but  is  typically 
of  order  of  one  hundred. 

In  Section  1.2.1,  necessary  notation  is  introduced  regarding  the  volume  and  projection  geometry  of 
polychromatic  cone  beam  diffraction  contrast  tomography  (henceforth  referred  to  as  LabDCT).  Section 
Section  1.2.2  and  Section  1.2.3  subsequently  introduce  the  back-  and  forward  projection  respectively. 
Finally,  in  Section  1.2.4,  we  discuss  some  implementation  details  that  were  required  to  optimise  to 
projector  performance. 

For  the  sake  of  readability,  we  refrain  from  delving  too  deep  into  the  mathematics  of  the  projection 
models.  Instead,  we  refer  to  Appendix  A  where  more  in-depth  derivations  are  provided  for  many  of  the 
formulas  presented  in  this  section. 

1.2.1  Notation 

Any  tomographic  projection  model  is  defined  by  both  a  volume  geometry,  defining  the  reconstruction 
volume,  and  a  projection  geometry,  defining  the  trajectory  and  orientation  of  the  X-ray  source  and 
detector1  with  respect  to  the  volume.  In  Fig.  1.3,  a  schematic  overview  is  given  of  a  lab  source  CT 
system  for  diffraction  tomography. 

For  the  volume  geometry,  define  n  as  the  total  number  of  voxels  and  Av  £  N  as  the  voxel  length  per 
direction.  For  simplicity’s  sake,  we  assume  isotropic  voxels. 

1  In  the  remainder  of  the  presentation  we  will  refer  to  individual  detector  pixels  as  detectors  and  we  will  use  the  term 
detector  system  when  referring  to  the  full  array  (2048x2048  elements)  of  detectors. 
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1.2.  CONE  BEAM  DCT  PROJECTOR 


Figure  1.2:  Schematic  showing  the  laboratory  DCT  setup.  The  technique,  originally  developed  for  the  case  of  a 
monochromatic,  parallel  synchrotron  beam  can  be  adapted  to  conventional  X-ray  laboratory  tomography  scanners 
with  minimal  modifications:  a  set  of  adjustable  slits  is  inserted  between  source  and  sample. 


For  the  projection  geometry,  we  define  m  as  the  total  number  of  detectors  per  projection  and  l  as 
the  total  number  of  projections.  The  total  number  of  measurements  is  then  fh  =  ml.  Let  At  £  R  denote 
the  width  and  height  of  the  square  detectors.  Each  projection  is  defined  by  a  set  of  five  vectors:  let 
5eR3  denote  the  position  of  the  source  point,  D0  £  R3  the  position  of  the  centre  of  the  detector  plane, 
U  £  R3  the  unit  vector  that  defines  the  x-dimension  of  a  detector  (||[/||  =  At),  V  £  R3  the  unit  vector 
that  defines  the  y-dimension  of  a  detector  (||Vj|  =  At),  and  G  £  R3  the  grain  lattice  plane  normal  giving 
rise  to  the  diffraction  spot.  Note  that  any  detector  point  can  be  expressed  by  D  =  D0  +  t\U  +  t2V  with 
ti,t2  £  N. 


1.2.2  Backprojection 

For  the  backprojection  operation  in  LabDCT,  only  a  voxel-driven  projection  approach  can  result  in  the 
fastest  implementations.  This  means  that  for  each  voxel  centred  at  point  P ,  and  for  each  projection 
image  defined  by  S  and  D0,  we  must  simulate  the  diffraction  of  the  incoming  ray  and  compute  the 
intersection  D  of  the  diffracted  ray  with  the  detector  plane  (see  Fig.  1.4).  Let  the  inverse  incoming  ray 
be  defined  by  the  vector  Ki„  =  P  —  S.  Given  the  law  of  ray  reflection,  the  outgoing  ray,  Kout  can  then 
found  by 

Kout  =  2G(G  ■  —Kin)  +  Kin.  (H) 

The  equivalence  of  Eqn.  1.1  to  Bragg’s  law  is  given  in  Appendix  A. 
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1.2.  CONE  BEAM  DCT  PROJECTOR 


Figure  1.4: 


Scheme  of  diffraction  in  lab  source  ct,  as  used  for  backprojection.  Simplified  to  2d  fan  beam. 


The  detector  D  can  then  be  found  by  computing  the  intersection  between  the  outgoing  ray  and  the 
detector  plane: 

P  +  sKout  =  D  =  Dq  +  t\U  +  t2V,  (l.2) 

where  s,  U  ,  t2  £  R  can  be  found  by  solving 

/ Kout,x  —Ux  — 14,\  /  s  \  / Dox  —  Px\ 

KOUt,y  -Uy  -Vy  k  =  \D0v  -  Py  ■  (H) 

\Kout,z  -Uz  -VZJ  \t2J  \D0z-PzJ 

Detailed  information  on  how  to  efficiently  solve  Eqn.  1.3  can  be  found  in  appendix  A. 

1.2.3  Forward  projection 

Here,  a  forward  projection  operator  is  discussed  for  use  in  LabDCT.  For  computational  reasons,  forward 
projecting  requires  a  ray-driven  projector.  Given  a  certain  detector  point,  we  must  thus  figure  out  which 
voxels  in  the  volume  contribute  exactly  how  much  to  it.  For  conventional  absorption  CT,  this  is  easy 
to  figure  out  as  all  contributing  points  lie  on  a  straight  line.  Rather,  for  diffraction  tomography,  these 
points  (i.e.  those  in  which  the  incoming  ray  diffracts  exactly  into  the  direction  of  the  detector)  are  on 
a  curved  line  (henceforth  referred  to  as  the  diffraction  curve),  requiring  more  complicated  projection 
equations. 

In  Section  1.2. 3.1,  we  will  introduce  formulas  that  define  the  diffraction  curve  for  the  special  case 
where  the  source  and  a  given  detector  are  on  the  cc-axis  of  a  2-dimensional  coordinate  system,  and  are 
located  such  that  their  distance  to  the  origin  is  equal.  Given  these  severe  restrictions,  we  can  make  use 
of  the  property  of  ellipses  saying  that  any  ray  originating  in  one  focal  point,  reflecting  on  the  edge  of  the 
ellipse,  will  end  up  in  the  other  focal  point. 

Subsequently,  in  Section  1.2. 3. 2,  a  transformation  formula  is  introduced  that  allows  the  special  ellipse- 
based  diffraction  curve  to  be  used  in  a  general  3D  LabDCT  cone  beam  setting. 
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1.2.  CONE  BEAM  DCT  PROJECTOR 


Figure  1.5:  A  subset  of  the  ellipses  with  their  centre  in  the  origin  and  with  fixed  focal  points  /  =  ±200.  The  thick 
red  line  connects  all  points  on  these  ellipses  where  the  corresponding  tangent  line  is  perpendicular  to  the  plane  normal 
(the  pink  line). 


1.2. 3.1  Ellipse  diffraction  model 

Consider  the  set  of  ellipses  with  its  centre  in  the  origin,  their  axes  aligned  with  those  of  the  Cartesian 
grid  and  with  fixed  focal  points  ±/  with  /2  =  a2  —  b 2  and  0  <  b  <  a. 


2  2 

Zi  ,  y_  = 

a2  b2 


(S) 


Given  the  in  depth  derivation  given  in  appendix  A,  we  can  see  that  all  points  x,  y  that  lie  on  such  an 
ellipse  and  whose  tangent  line  at  that  point  of  the  ellipse  is  perpendicular  to  the  diffraction  normal  G , 
can  be  described  as  a  function  of  a: 


a2\Gx\  _  sgn  (Gx)b2Gy 

x  v^7 

The  diffraction  curve  can  also  be  described  as  a  function  of  x: 

Gy  2  f2GXGy 

x±  Jx2  -  4/2G2G2 


HD 


HD 


In  Fig.  1.5  an  example  set  of  ellipses  is  shown  that  conform  to  Eqn.  1.4  with  /  =  ±200.  It  also  shows 
all  points  x,y  of  the  diffraction  curve  (Eqn.  1.5)  for  a  £  R+. 


1.2. 3. 2  Cone  beam  geometry 

Now  that  we  know  the  diffraction  curve  for  the  case  where  the  source  and  detectors  are  located  on 
the  focal  points  of  ellipses  centred  on  the  origin,  we  introduce  a  transformation  such  that  each  source, 
detector  pair  in  volume  coordinates  is  convert  into  ellipse  coordinates  such  that  these  restrictions  apply. 

Given  are  a  source  point  S  and  detector  point  D.  The  focal  points,  ±/,  of  all  ellipses  is  half  the 
Euclidean  distance  between  source  and  detector: 

/  =  \\SDh  =  i_^/{Dx  _  Sx)2  +  {Dy  _  Sy)2  +  {Dz  -  Szy.  dz) 

With  /  known,  the  source  and  detector  points  can  be  expressed  in  ellipse  coordinates: 

s'  =  (-f,  o)  £>'  =  (/,  0).  dD 
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1.2.  CONE  BEAM  DCT  PROJECTOR 


Figure  1.6:  Two-dimensional  simplification  of  the  transformation  between  the  a  volume  coordinates  (full  lines)  and 
ellipse  coordinates  (dotted  lines). 


The  plane  normal  in  ellipse  coordinates,  G' ,  can  be  described  by  the  angle  9  between  G  and  the  vector,  B, 
from  the  halfway  point  between  the  source  and  the  detector,  to  the  detector  (i.e.  B  =  D  — 

From  the  definition  of  the  dot  product,  we  learn  that 

„  BG  BG 

cos^MM~T 

Given  sin2  9  +  cos 2  9  =  1,  sinf?  can  be  also  be  computed  without  performing  a  costly  trigonometric 
function.  The  plane  normal  in  ellipse  coordinates  is  thus: 

G'  =  \/l  —  g ,  with  g  =  —  ( BXGX  +  ByGy  +  BZGZ) .  (1-10  j 

If  it  is  assumed  that  G'  is  not  parallel  to  D'  (which  means  the  ray  will  never  hit  the  detector  anyway), 
the  vectors  G '  and  D '  provide  a  basis  for  the  ellipse  coordinate  system  and  each  point  P'  can  be  described 
by  a  linear  combination  of  these  vectors: 


P’  =  olD’  +  /3  G1 


P' 


P 


af  +  pG'x,  P'y  =  f3G'y 
P y_  _  PL-  PG’X 

G'y  f 


Furthermore,  note  that  the  transformation  of  G'  and  D'  into  volume  coordinates  (i.e.  G  and  B)  provides 
a  basis  for  all  diffraction  points  in  volume  coordinates  with  the  origin  at  O  =  SL2D .  Because  the 
transformation  from  ellipse  to  volume  coordinates  involves  translation  and  rotation,  but  no  scaling,  the 
linear  combination  of  these  base  vectors  in  the  ellipse  coordinate  system  is  equal  to  the  corresponding 
set  of  base  vectors  in  the  volume  coordinate  system,  i.e. 

P  =  0  +  aB  +  /3G.  (Tl2j 


This  means  that  we  now  have  the  means  to  transform  the  LabDCT  projection  of  any  source-detector 
pair  into  the  simplified  geometry  of  Section  1.2. 3.1  (with  Eqn.  1.10),  where  the  diffraction  curve  can 
be  easily  computed.  With  Eqn.  1.12,  this  curve  can  be  reconverted  into  the  actual  volume  coordinates, 
leading  to  a  forward  projection. 
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1.2.  CONE  BEAM  DCT  PROJECTOR 


1.2. 3. 3  Sampling  based  projector 


Figure  1.7:  Scheme  depicting  a  two-dimensional  simplification  of  volume  sampling  along  a  certain  diffraction  curve. 

Now  that  formulas  for  the  diffraction  curve  are  available,  we  can  use  them  to  create  a  ray-driven 
forward  projection  strategy.  We  have  found  that  from  a  computational  point  of  view,  the  best  projection 
method  is  to  sample  the  volume  along  the  diffraction  curve  Eqn.  1.6  at  a  certain  interval  (see  Fig.  1.7). 
With  such  methods,  the  sampling  interval  should  be  as  regular  as  possible.  In  the  case  of  the  diffraction 
curve  however,  the  derivative  of  the  curve  is  not  constant  and  if  it  were  to  be  sampled  with  a  fixed  A.t 
or  a  fixed  Ay,  it  is  likely  to  be  severely  undersampled  when  its  slope  is  very  steep. 

Instead,  note  that  for  each  value  of  x  (with  f  <  x  <  —  \2fGxGy\),  a  unique  value  of  a  can  be 
associated,  and  that  a  increases  as  x  increases  (This  is  only  true  for  x  <  0.  However,  because  we  can 
assume  that  the  reconstruction  volume  is  much  closer  to  the  source  than  it  is  to  the  detector,  the  case 
where  x  >  0  can  be  ignored).  One  could  thus  also  sample  the  curve  by  slowly  increasing  Aa,  which  offers 
a  much  more  even  sampling  strategy.  This  has  the  added  benefit  that  the  projection  can  be  computed 
using  Eqn.  1.5,  which  is  substantially  more  computationally  efficient  than  Eqn.  1.6.  2 

Given  0  <  ao  <  /  a  certain  starting  value,  compute  its  corresponding  point,  Pq,  on  the  diffraction 
curve.  Next,  add  Aa  to  ao  and  compute  its  corresponding  point  Pi.  Before  Pi  is  "looked”  up  in  the 
volume  image  and  added  to  the  detector  value,  its  weight  has  to  be  computed.  This  weight  should 
be  the  inverse  of  the  curve  length  between  P0  and  Pi.  Because  it  is  not  possible  to  compute  this 
efficiently,  one  simplification  can  be  to  use  the  Euclidean  distance  between  the  two  points.  If  Aa  is 
chosen  sufficiently  small,  this  simplification  should  not  result  in  large  errors.  Keep  adding  Aa  to  a  and 
repeat  the  computations  until  the  diffraction  curve  does  not  lie  in  the  volume  any  more. 

As  an  additional  improvement,  one  can  adaptively  update  Aa  by  dividing  it  by  the  distance  between 
the  last  two  points  (i.e.  the  weight  associated  to  the  last  point).  This  will  ensure  that  the  distance  from 
this  point  to  the  next  will  be  approximately  equal  to  the  distance  between  this  point  and  the  previous, 
obviously  resulting  in  an  a  much  smoother  sampling  of  the  curve. 

In  Fig.  1.8,  the  LabDCT  forward  projector  is  summarized  in  pseudo  code. 

1.2. 3. 4  Some  optimizations 

We  have  done  research  into  the  further  improvement  of  the  performance  of  the  forward  projection  opera¬ 
tion.  We  noticed  that  in  the  typical  setup,  the  reconstruction  volume  lies  much  closer  to  the  source  than 
to  the  detector  and  is  very  small  relative  to  the  source  detector  distance.  Furthermore,  for  grains  to  be 
visible  on  the  detector  plane,  the  plane  normal  of  the  projection  is  typically  close  to  perpendicular  to  the 
incoming  rays.  This  means  that  the  curvature  of  the  portion  of  the  diffraction  curve  that  falls  inside  the 
region  of  interest,  is  typically  very  low.  We  have  thus  found  that  the  curve  can  be  approximated  with  a 
straight  line  without  a  substantial  loss  in  accuracy. 

2Note  that  in  Eqn.  1.6  the  computational  burden  appears  to  lie  mostly  with  the  reciprocal  square  root.  Fortunately, 
however,  on  nVidia  GPU’s  the  reciprocal  square  root  is  implemented  as  a  hardware  instruction,  and  thus  requires  only  one 
clock  cycle  to  compute. 
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1.3.  EXPERIMENT 


a  :=  ao  such  that  Po  lies  inside  the  volume 

PL  ■=  ~c?\Gx\l -  pG%-,  P'y  =  -sgn(Gx)b2\Gv\/^ai-  pGf, 

P-=n  y/G'yi  <*:=(PL-pG'x)/f-, 

P0:=O  +  cx.B  +  f3G', 
a  :=  a  +  A  a; 

while  Po  is  inside  the  volume 

P[x  ■=  —a2\Gx\/ yj  a?  —  pG^;  P[y  =  -sgn(Gx)P\Gy\/^/^  -  pG$; 
P  :=  PL/G'y,  oc~(PL-pG'x)/f-, 

Pi  :=  O  +  cxB  +  pG\ 
w  :=  ||Pi  -  -P0II2; 

v  :=  lookup  value  of  the  volume  at  point  Pi; 

projection  :=  projection  +  wv: 

Po  :=Pi; 

A  a  :=  Aa/w;  (optional) 
a  :=  a  +  A  a; 

end 


Figure  1.8:  Pseudo  code  for  sampling-based  LabDCT  forward  projector. 


Assume  that  the  grain  falls  entirely  inside  a  sphere  centred  in  point  C  G  R2  with  radius  r  G  R.  One 
can  then  compute  the  intersection  of  this  sphere  with  the  diffraction  curve  Eqn.  1.5.  One  can  then  draw 
a  straight  line  between  these  two  points,  sampling  them  at  a  fixed  interval. 

1.2.4  Implementation  details 

In  this  work,  multiple  implementations  of  the  previously  discussed  projector  were  created.  Initial  proto¬ 
typing  development  was  done  fully  in  MALTAB.  Later  on,  a  C++  implementation  was  created  where  the 
computational  efficiency  was  optimized.  In  this,  the  inherent  parallelism  of  the  projectors  was  exploited 
with  the  OpenMP  framework.  To  even  further  decrease  the  computation  time,  also  an  NVIDIA  CUDA 
implementation  was  created.  In  it,  parallelism  is  exploited  even  more  ad  GPU  specific  techniques  such 
as  texture  memory  are  used  to  increase  the  speed  of  bicubic  interpolation  and  random  memory  lookups. 

Both  the  C++  and  the  CUDA  versions  are  part  of  the  ASTRA  framework.  For  now,  access  to 
this  code  is  restricted  to  the  partners  of  this  project.  There  is,  however,  the  option  to  include  these 
projectors  in  a  future  open  source  release  of  the  ASTRA  toolbox.  At  the  time  of  writing,  the  code  is 
already  deployed  on  a  computation  cluster  at  the  ESRF,  Grenoble,  France. 

1.3  Experiment 

For  validation  of  the  projector  introduced  in  Section  1.2,  we  consider  a  experimental  dataset  acquired  on 
the  laboratory  X-ray  tomography  setup  at  MATEIS,  INSA  Lyon.  The  instrument  (manufacturer:  RX 
Solutions,  France)  is  equipped  with  a  Hamamatsu  microfocus  X-ray  source.  The  source  was  operated  at 
60  keV  acceleration  voltage  and  a  5/irn  thick  tungsten  transmission  target.  The  images  were  recorded 
on  a  CCD  detector  (ESRF  Frelon  camera)  featuring  a  2048  x  2048  pixel  array  and  a  dynamic  range  of 
14  bits.  The  backilluminated  CCD  sensor  is  fibre-optically  coupled  to  a  100  m  thick  fluorescent  screen 
(Gad Ox).  The  pixel  size  of  the  detector  system  was  50  microns.  The  centre  of  rotation  of  the  tomographic 
rotation  stage  was  positioned  at  a  distance  of  36  mm  from  the  source  and  the  detector  was  positioned  at 
a  distance  of  220  mm  from  the  source.  This  resulted  in  an  effective  pixel  size  of  about  8  microns,  taking 
into  account  the  cone  beam  magnification.  A  set  of  adjustable  slits  were  inserted  between  the  sample 
and  the  source  in  order  to  limit  direct  X-ray  exposure  to  the  central  part  of  the  detector  field  of  view 
(figure  1.2).  Two  scans  were  recorded  from  the  same  sample  —  a  DCT  scan  optimized  for  detection  of 
the  weak  diffraction  signals  and  a  conventional  absorption  scan.  For  the  DCT  scan,  a  1  mm  thick  copper 
sheet  was  placed  as  an  absorber  over  the  area  of  the  detector  exposed  to  the  direct  beam,  in  order  to 
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1.3.  EXPERIMENT 


Figure  1.9:  Slice  of  the  absorption  reconstruction  of  the  AlCu  dataset.  Highlighted  is  the  grain  that  is  used  in  the 
subsequent  experiments. 


avoid  saturation  of  the  detector  system. 

The  sample  was  made  from  a  binary  A1  Cu  alloy  (5  weight  percent  copper)  which  was  thermo- 
mechanically  processed  in  a  way  to  produce  large  A1  grains  (  500  microns)  partly  surrounded  by  a 
copper  rich  Eutectic  phase.  A  slice  of  the  absorption  reconstruction  of  this  dataset  is  shown  in  Fig.  1.9. 

One  reason  we  chose  this  dataset  is  that  the  individual  grains  can  be  more  or  less  seen  in  the  absorption 
reconstruction,  allowing  us  to  create  a  manual  three-dimensional  segmentation  of  the  grain  shapes.  Note 
that  while  it  is  tempting  to  refer  to  such  a  segmentation  as  a  ’’ground  truth” ,  the  fact  that  it  was  created 
by  manually  delineating  the  borders  of  the  grain  on  a  slice-by-slice  basis,  means  that  it  should  only  be 
used  as  a  rough  guideline  of  the  exact  grain  shape.  A  more  accurate  segmentation  may  be  obtained  using 
advanced  morphological  image  segmentation  techniques  (watershed  algorithm). 

Following  the  procedure  outlined  in  [4],  a  number  of  grains  were  identified  (indexed)  from  the  LabDCT 
dataset.  In  what  follows,  we  describe  two  experiments  that  we  did  with  the  data  from  one  of  those  grains 
which  could  be  identified  both  in  the  absorption  reconstruction  (i.e.  a  grain  delineated  by  the  Eutectic 
phase)  and  indexed  during  processing  of  the  diffraction  data. 

First,  in  Section  1.3.1,  we  show  a  reconstruction  that  we  did  of  a  phantom  volume,  based  on  the 
manually  segmented  grain  volume.  Secondly,  in  Section  1.3.2  we  attempt  to  reconstruct  the  same  grain, 
with  the  actual  experiment  data. 

1.3.1  Phantom  experiment 

Firstly,  we  consider  the  manual  segmentation  we  created  of  the  grain  highlighted  in  Fig.  1.9,  a  slice  of 
which  is  shown  in  Fig.  1.10a.  Of  this  volume,  we  simulated  the  diffraction  spots  with  the  same  projection 
geometry  (i.e.  the  same  directions)  as  the  experimentally  recorded  dataset.  With  this  projection  data, 
we  computed  a  reconstruction,  shown  in  Fig.  1.10b.  As  reconstruction  algorithm,  we  used  50  iterations  of 
SIRT  with  minimum  constraints  (i.e.  after  each  SIRT  iterations,  all  negative  values  in  the  reconstructed 
are  put  to  zero). 

While  there  are  some  clearly  visible  streaks  in  the  reconstructed  image,  the  general  outline  of  the 
grain  is  reconstructed  very  well. 

1.3.2  Real  data  experiment 

Secondly,  we  created  a  SIRT  reconstruction  of  the  experimental  data.  This  time,  we  opted  for  only 
10  iterations,  as  we  noticed  a  decreasing  reconstruction  quality  when,  due  to  noisy  projection  data, 
doing  more  iterations.  A  slice  of  the  reconstruction  is  shown  in  Fig.  1.11b,  where  it  is  also  compared 
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1.3.  EXPERIMENT 


(a)  manual  segmentation  (b)  reconstruction 


Figure  1.10:  (a)  Slice  of  the  manually  segmented  grain  (Fig.  1.9)  (b)  Reconstruction  of  (a)  computed  from  162 
simulated  diffraction  spots  with  the  same  geometry  as  the  dataset  of  Fig.  1.9. 


(a)  approximation  [4]  (b)  LabDCT  (Section  1.2)  (c)  approximation  [4]  (d)  LabDCT  (Section  1.2) 

Figure  1.11:  (a)  Reconstruction  slice  with  an  approximated  projector  [4],  (b)  Reconstruction  slice  with  the  projector 
introduced  in  Section  1.2.  (c,d)  Volume  rendering  of  the  full  volumes  of  (a)  and  (b)  respectively. 


(Fig.  1.11a)  to  the  reconstruction  obtained  with  the  previously  published  approach  in  which  regular  cone 
beam  projectors  were  used  and  the  diffraction  geometry  was  approximated  by  some  afffine  transformation 
of  the  projection  images  [4]. 

Clearly,  the  method  presented  in  this  work  much  better  approximates  the  shape  of  the  grain  as  it  is 
visual  in  the  absorption  reconstruction  (Note  that  the  shown  slices  in  Fig.  1.11a,  Fig.  1.11b  and  Fig.  1.9 
are  not  exactly  the  same). 

We  also  compared  the  computational  efficiency  of  the  different  implementations  of  the  projectors. 
In  Fig.  1.12,  timings  are  shown  for  the  forward  and  backprojection  and  for  a  10  iteration  SIRT  recon¬ 
struction.  The  computations  were  performed  on  an  Intel  Quad  Core  i7-2760  CPU  running  at  2.40GHz 
and  on  an  NVIDIA  Quadro  2000M  GPU  (note  that  this  is  only  a  laptop  GPU,  with  a  modern  high 
end  GPU  we  expect  a  further  speedup  by  a  factor  of  2).  An  enormous  computational  improvement  was 
obtained  by  implementing  the  projector  in  C++  instead  of  in  MATLAB  code.  While  this  was  indeed 
expected,  it  must  be  noted,  however,  that  the  MATLAB  implementation  is  not  fully  optimized  for  speed 
and  the  comparison  is  thus  unfair.  Overall,  however,  it  is  clear  that  for  the  best  performance,  GPU 
computing  is  the  way  to  go.  With  features  such  as  fast  texture  interpolation  and  memory  lookup,  the 
GPU  architecture  is  excellently  suited  for  projectors  such  as  ours,  and  vice  versa. 
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implementation 

FP 

BP 

SIRT  (10  iters) 

MATLAB 

896.40s 

2128.12s 

32514.11s 

C++ 

5.89s 

28.10s 

382.64s 

CUDA 

0.45s 

0.69s 

13.78s 

Figure  1.12:  Projection  and  reconstruction  timings  for  the  different  implementations. 

1.4  Outlook  and  conclusions 


In  this  work  we  have  developed  an  efficient  projector  for  laboratory  based,  polychromatic  cone  beam 
diffraction  imaging.  Using  this  projector,  the  three-dimensional  shape  of  grains  with  negligible  intra- 
granular  orientation  gradient  can  be  reconstructed  using  iterative,  algebraic  reconstruction  techniques 
(e.g.  SIRT).  The  reconstructions  show  a  clear  improvement  compared  to  the  previously  used  processing 
route,  based  on  a  regular  cone  beam  projector  and  affine  transformation  of  the  projection  images. 

Further  validation  and  evaluation  of  the  new  projector  will  become  possible  after  sucessfull  indexation 
of  the  majority  of  the  grains  in  the  illuminated  sample  volume  -  currently  only  about  40  percent  of  the 
sample  volume  has  been  indexed.  Where  this  document  only  provides  a  visual  confirmation  of  the 
accuracy  (partly  delineation  of  the  grain  boundaries  by  the  Eutectic  phase),  quantiative  evaluation  will 
be  possible  after  having  scanned  the  same  sample  on  a  synchrotron  beamline  using  the  well  established 
monochromatic,  parallel  beam  approach. 

We  expect  that  further  improvement  of  the  reconstruction  quality  can  be  obtained  from  more  ac¬ 
curate  estimation  of  the  projection  geometry.  Currently,  the  plane  normals,  defining  the  direction  of 
the  diffracted  beam  are  derived  from  the  orientation  matrix,  determined  during  the  indexing  procedure. 
Due  to  incertainties  in  the  global  setup  parameters,  this  initial  estimate  of  the  orientation  is  not  100 
percent  accurate,  which  results  in  blurry  reconstructions.  This  problem  is  closely  related  to  projection 
misalignment,  common  to  other  applications  such  as  electron  tomography. 

The  logic  continuation  of  the  current  work,  aiming  at  further  improvement  of  the  grain  shape  recon¬ 
structions  is  along  the  following  lines: 

Grain  refinement  The  initial  estimate  of  grain  center  of  mass  position  and  orientation  as  provided 
by  the  indexing  routine  has  to  be  refined.  However,  the  refinement  of  grain  parameters  can  only 
be  successful  if  the  global  setup  parameters  like  detector  position  and  tilt  of  the  rotation  axis  are 
refined  at  the  same  time.  Therefore,  a  multiparameter  fitting  procedure,  based  on  the  full  set  of 
diffraction  spot  metadata  is  required  (estimated  time:  1-2  months  for  implementation  and  testing). 
Alternatively,  grain  reconstructions  could  be  improved  on  a  grain  by  grain  basis.  Recently,  at  the 
Vision  Lab  (Antwerp,  Belgium),  a  state-of-the-art  alignment  correction  method  for  tomography 
was  developed  (not  yet  published).  We  are  confident  that  the  underlying  principles  of  this  method 
can  also  be  used  to  improve  the  plane  normal  estimations,  which  is  likely  to  lead  to  sharper  recon¬ 
structions  of  individual  grains,  albeit  not  solving  the  underlying  problem  of  a  slightly  inaccurate 
initial  position  and  orientation  estimate  of  the  grain. 

Intragranular  orientation  gradients  One  of  the  principal  limitations  of  polychromatic  cone  beam 
DCT  is  its  restriction  to  materials  exhibiting  negligible  levels  of  intragranular  orientation  spread. 
In  the  current  work  we  have  assumed  that  the  direction  of  the  diffracting  plane  normal  is  constant 
throughout  the  entire  grain  volume.  In  practise,  this  condition  is  rarely  fulfilled  and  already 
small  deviations  of  order  of  a  few  tenths  of  a  degree  of  intragranular  orientation  spread  can  give 
rise  to  significant  distortion  of  the  diffraction  spots  and  will  lead  to  rapid  degradation  of  the 
accurary  of  grain  shape  reconstructions.  It  can  be  expected  that  a  6-dimensional  reconstruction 
framework,  as  currently  developed  for  the  case  of  monochromatic  beam  DCT  [5],  will  be  able 
to  improve  the  quality  LabDCT  reconstructions,  as  well.  Since  this  method  is  computationally 
expensive,  the  fast  GPU  implementation  developed  in  the  current  work  may  be  considered  as  a 
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prerequiste  and  first  step  in  this  direction.  However,  the  polychronraticity  of  LabDCT  implies  that 
the  additional  information  encoded  into  the  u>  (rotation  angle)  spead  of  the  individual  diffraction 
spots,  observed  in  the  monocliroamtic  beam  variant,  is  lost.  As  a  first  step,  the  stability  and 
applicability  of  such  a  6-D  reconstruction  framework  has  to  be  tested  on  synthetic  data  created 
from  a  phantom  grain  with  known  orientation  distribution.  If  successfull  the  application  to  real  data 
and  validation  with  cross-sectional  EBSD  and  /  or  synchtrotron  measurements  can  be  envisaged. 
The  time  required  for  adapting  the  6-dimensional  reconstruction  framework  to  LabDCT  (creation 
of  test  cases,  optimization  of  the  algorithm  and  model  ,  validation)  is  estimated  to  be  of  order  of 
6  months. 

Source  spectrum  In  this  work,  we  have  assumed  that  the  polychromatic  X-ray  source  has  a  constant 
energy  spectrum.  In  practice  this  may  not  be  the  case  and  strong  peaks  (characteristic  emission 
lines  of  the  X-ray  target  material)  can  be  present  in  the  spectrum.  This  means  that  within  a  single 
diffraction  spot,  which  typically  covers  a  small  range  of  different  energies,  the  contribution  from 
different  grain  sub-volumes  may  be  highly  non-uniform.  Assuming  the  source  spectrum  is  known 
a  priori,  one  may  attempt  to  extend  the  current  projection  model  such  that  this  non-constant 
spectrum  is  taken  into  account.  While  this  is  likely  to  decrease  the  computational  efficiency  of 
the  projection,  it  may  be  beneficial  to  its  accuracy.  As  first  step,  we  intend  to  perform  software 
simulations  in  order  to  evaluate  the  impact  of  the  source  spectrum  on  the  reconstruction  quality.  If 
these  simulations  indicate  a  strong  influence,  we  will  proceed  with  the  implementation  of  projectors 
and  algorithms,  taking  this  effect  into  account. 

To  summarize,  we  have  implemented  and  validated  a  projector  for  LabDCT  which  enables  us  to  per¬ 
form  accurate  3D  grain  shape  reconstructions  in  materials  exhibiting  negligible  intra-granular  orientation 
spread.  With  further  improvement  of  polycrystal  indexing  algorithms  and  with  extension  to  a  6D  re¬ 
construction  framework,  we  anticipate  LabDCT  to  become  a  true  alternative  to  established  synchrotron 
techniques.  The  combination  of  absorption  (detection  of  damage,  second  phases,  porosity,  plastic  strain 
via  3D  digital  volume  correlation  techniques)  and  diffraction  imaging  (crystal  shape  and  orientation, 
elastic  strain)  on  the  same  instrument  will  provide  unique  possibilities  for  (interrupted)  in-situ  charac¬ 
terization  of  dynamic  processes  in  the  bulk  of  polycrystalline  materials  (e.g.  grain  coarsening  during 
heat  treatment,  plastic  deformation  and  damage  mechanisms,  others...) 
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2 

Appendix 


2.1  Backprojection 

2.1.1  Link  to  Bragg’s  law 


Figure  2.1:  Diffraction  direction  according  to  Bragg’s  law. 


In  DCT,  the  direction  of  the  diffraction  is  defined  by  the  Laue  equation,  which  is  equivalent  to 
Braggs’s  Law: 

G  =  Kout  -  Kin,  (2T) 

where  ||A'out||  =  ||Ajn||,  and  ||AT ||  =  1  (see  Fig.  2.1).  Here,  we  derive  the  equivalence  to  Eqn.  1.1.  Given 
the  incoming  angle  a ,  the  norm  of  Kout  and  Kin  can  be  easily  found: 


We  then  get: 


1  =  ||G||  =  2||/Lm||  cos ct 


\Kin\\  = 


1 


2  sin  6 


2\\Kin\\  sin 9, 


Kout  =  2 G(G  ■  Kin )  +  Kin 

=  2G{\\G\\\\Kin\\sme)+Km 

=  2G(n  1  Q  Sin6>)  + 

2  sm^ 

=  G+  Kin  @ 
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2.2.  FORWARD  PROJECTION,  ELLIPSE  MODEL 


2.1.2  Solving  Eqn  1.3 

Given  the  equation  Dq  +  t\U  +  £2^  =  P  +  sKout,  here  we  discuss  how  to  efficiently  find  s,tl,t2  £  R. 
First,  find  s  by  computing  the  intersection  between  the  line  P  +  sKout  and  the  detector  plane: 

=  (A)  -P)-N 
S  Kout  ■  N  ’ 

where  N  is  the  normal  to  the  detector  plane,  i.e.  N  =  U  x  V .  Note  that  N  can  be  precomputed  for 
computational  efficiency. 

The  value  for  1 1  can  subsequently  be  computed,  depending  on  the  value  for  Ux: 

Vz(Px-Dox-\-sKouttX')  —  Vx(.IDz—Da2-\-8KouttZ)  •  y  jj  /  r> 

UXVZ-U,VX  11  Ux  7=  u 

VZ{Py  —  DQy-\-SKout  y)—Vy{PZ—DQz+SKout^z)  \(  TT  '  O 
UyVZ-UZVy  ii  Uy  ^  U 


dD 


Note  that  we  don’t  consider  the  case  where  Ux  =  Uy  =  0.  However,  because  U  represents  the  unit  vector 
for  a  column  on  the  detector  plane,  we  can  assume  that  at  least  one  of  Ux  and  Uy  is  non-zero.  For  the 
same  reason  we  can  also  assume  that  Vz  ^  0,  leading  to  the  following  equation  for  t2: 

,  _  Pz  -  Dqz  +  s  *  Kout,z  —  tl*Uz 


2.2  Forward  Projection,  Ellipse  model 


Consider  the  set  of  ellipses  with  its  center  in  the  origin,  their  axes  aligned  with  those  of  the  Cartesian 
grid  and  with  fixed  focal  points  ±/  with  /2  =  a2  —  b2  and  0  <  b  <  a. 


b2 


=  1. 


Given  a  certain  point  Xo,  yo  on  such  an  ellipse,  the  slope  of  the  tangent  line  of  the  ellipse  at  that  point 
can  be  found  by  implicit  differentiation  along  x: 


b2x2  +  a2y2  =  a2b2 
’I,  (b2x2  +  a2  +  y{ x)2)  =  0 
2  b2x  +  2  a2y^~  =  0. 


The  slope  of  the  tangent  line,  mt,  is  thus 

dy  b2  x 

mt  =  -1-  = - 2  * 

ax  az  y 

We  are  now  looking  for  points  Xo,yo  where  the  tangent  line  is  perpendicular  to  the  diffraction  normal, 
whose  slope  is  defined  by  mg  = 


TOtTOg  =  -1 
b2  X0  Gy  _ 

a 2  y0  Gx 
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2.2.  FORWARD  PROJECTION,  ELLIPSE  MODEL 


This  gives  rise  to  a  line  through  the  origin,  that  intersects  the  ellipse  exactly  at  the  two  points  where  its 
tangent  curve  is  perpendicular  to  the  normal  vector: 


a2Gxy  =  b2Gyx 


ti^Gy  b2 
a2  (L' 


Now,  given  a  certain  Xq  ,  we  wish  to  know  yo  for  which  the  previous  equations  hold.  This  can  be  done 
by  computing  the  intersection  between  Eqn.  2.9  and  the  ellipse  Eqn.  2.7,  which  is1: 


x0  =  ± 


ab 


=  ± 


a2^#+62 


ab  b2 


b2m 2 


=  ± 


2  _  /!"*!  ,  x 

mg  a2  +  1 


=  ± 


1  G2yf 2 
G2  G2a2 


=  ± 


a2\Gx 


a2b2\Gx\Gy  sgn  (Gx)b2Gv 

yo  =  ±  7  ,  =z?m9  =  ± - /  =  ±  ,  ■ 

Gxa\a?-pG*  Ja2-f2G2y 


a2  b~xrL  +  b2  a 


}/* 2~f2Gl 

o 

(HD 


Given  that  the  volume  will  always  lie  much  closer  to  the  source  than  to  the  detector,  we  are  only  interested 
in  points  for  which  Xq  <  0.  Thus, 


1 


7  = 


- =,  x0  =  a2\Gx\i,  y0  =  sgn(Gx)Gy(f2  -  a2)'). 

\/«2  -  pgi 


(HD 


With  xq  given  and  b  easily  computable  with  b2  =  a2  —  / 2,  Eqn.  2.12  can  be  solved  for  the  two  remaining 
unknowns  (yo  and  a).  From  Eqn.  2.10,  a 2  can  be  found: 


2  a‘ G2X 

X°  ~  a2-  f2G2 
0  =  G2xaA  -  x20a2  +  f2G2yx20 

Solving  for  a 2,  one  gets 

°2  =  2^|  (*°  *  xo\J xo  ~  ^G2G2f2^J  .  \2A3j 

With  a2  known  and  b2  =  a2  —  /2,  Eqn.  2.9  can  be  used  to  represent  the  diffraction  curve  yo  as  a  function 
of  Xq’- 


b2  f2 

VO  =  a2m9X°  =  m9X°  ~  ^2m9X0 

_  Gy  r _ f2GyX  Q _ 

G"  a°  m  (xo  ±  \JXq  —  4G2G2/2)  G, 
Gv  2  f2GxGy 

=  yr-XO - 7  =■ 

Xo  ±  ^x2  -  Af2G2G2 

It  should  be  noted  that  if 

x20  <  4/2G2G2  |x0|  <  \2fGxGy\, 

1  http:  /  /  mathworld.wolfram.com /Ellipse-Linelntersection.html 


(HD 


(HD 
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2.2.  FORWARD  PROJECTION,  ELLIPSE  MODEL 


a2  will  be  a  complex  number,  meaning  that  for  these  Xo  values,  no  ellipse  exists  where  the  tangent  slope 
at  xq  is  perpendicular  to  G.  This  property  can  be  used  to  define  limits  to  where  the  curve  should  be 
computed. 

Also  note  that  both  Eqn.  2.12  and  Eqn.  2.13  feature  <4-.  This  is  a  potential  cause  for  severe  floating 
point  errors  if  Gx  is  very  small.  In  that  case,  the  plane  normal  is  close  to  perpendicular  to  the  incoming 
ray  and  the  diffraction  will  be  infinitesimally  small. 
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